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SOLUTIONS OF DIFFERENTIAL EQUATIONS WITH REGULAR COEFFICIENTS 
BY THE METHODS OF RICHMOND AND RUNGE-KUTTA 


ABSTRACT 


Numerical solutions of the differential equation which 
describes the electric field within an inhomogeneous layer of 
permittivity, upon which a perpendicularly-polarized plane wave is 
incident, are considered. Richmond's method and the Runge-Kutta 
method are compared for linear and exponential profiles of 
permittivities. These two approximate solutions are also compared 
with the exact solutions. 
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INTRODUCTION 


The treatment of perpendicularly-polarized plane wave 
incidence on inhomogeneous dielectric or plasma layers has been 
investigated by many authors [1-3]. Their efforts have been 
directed toward solving the differential equations which describe 
the electromagnetic behavior of fields within these layers. 
Solutions for continuous inhomogeneous layers require, in general, 
numerical techniques. Exact solutions for linear and exponential 
profiles are available [4]. 

Many methods have been used to solve the differential 
equations describing the fields within an inhomogeneous layer. 
Dividing the layer into a number of smaller segments and 
approximating the inhomogeneous profile by constant permittivity 
values within each segment, the multilayer boundary value problem 
is then solved. Approximate integral solutions are also available 
as are WKB asymptotic solutions where the inhomogeneous profile 
(or permittivity) gradient is small [5,4]. For these solutions 
the computation time is too long or the model is too crude to 
yield accurate results. 

Practical solutions have been obtained using the Runge-Kutta 
method [6] and a step-by-step numerical method described by 
Richmond [7]. Both solutions have been programmed on a personal 
computer and the results compared. A comparison with exact 
solutions, linear and exponential, has also been made. 

The primary purpose of this paper is to explain and 
demonstrate the Richmond and Runge-Kutta methods for solving a 
differential equation in which its parameters may or may not be 
known analytically. For definiteness, therefore, the differential 
equation which represents the electric field within a plane 
inhomogeneous layer of dielectric or plasma, upon which a 
perpendicularly-polarized plane wave has been assumed incident, 
has been considered. An efficient solution of this differential 
equation for this ideal model is needed since the same 
differential equation exists for realistic problems such as 
aperture antennas which are bounded by inhomogeneous dielectric or 
plasma layers. 



INHOMOGENEOUS PLANE LAYER-PERPENDICULAR POLARIZATION 


A plane wave is incident on a plane inhomogeneous layer of 
dielectric or plasma as shown in figure 1. The incident electric 
field which is polarized in the x-direction propagates in free 
space (region I) at an angle 8 with the normal to the layer. The 
transmitted electric field which is also x-polarized exists into 
free space (region III) at an angle 8. The partial differential 
equation for the electric field in the free space regions and in 
the inhomogeneous layer is written as 

d 2 E (y,z) a 2 E (y,z) 2 

— 2 + — 2 ^ + k c r (z) E X (Y ' Z) = 0 (1) 

3y 3z 

where k is the free space wavenumber and e r (z) is the relative 

dielectric constant of the medium. Since the relative dielectric 
constant in free space (regions I and III) equals one, the 
solutions of equation (1) in these regions are readily found: 

E I (y,z) = E ± [ e jkzcos ® + R e'jkzcosaj e j*ysine (2) 

E* Ir (y,z) = E t ei kzc ° se e jkysine (3) 

where E^ is the intensity of the incident electric field, R is the 

reflection coefficient at the interface, and E fc is the transmitted 

electric field intensity. However, in region II, since e r (z) can 

be an arbitrary function of z, the solution of equation (1) must 
be determined numerically, except for some specialized profiles 
such as linear and exponential. 

Before attempting the solution of equation (1) in region II, 
it is reduced to an ordinary differential equation. This is 
readily accomplished by the method of separation of variables and 
the continuity of tangential fields throughout the media. The 
continuity requirement forces the y-variation in the fields to be 
the same as given in equations (2) and (3) . The differential 
equation describing the z-variation of the electric field is 
written as 

d 2 E 11 (z) 

+ k 2 ( e (z) - sin 2 e ) E^ 1 (z) = 0 (4) 

dz z r x 
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Richmond ' s Method 


One method of solving equation (4) for region II is reported 
in a paper by Richmond [7], In his method a step-by-step 
numerical procedure is used to find the electric field within the 
layer at successive points. To begin the solution, the field at 
z=0 and z=h, which is a small distance within the layer, must be 
determined first. The field at z=0 is assumed to be unity and the 
field at z=h is found from the Maclaurin series expansion at 
z=0; i.e. , 


E^(h) 


where 


dE x (z) 


= E ( 0 ) + -3- 

x v ' dz 


V° ) = 1 


h + _L_ 

d E x (z) 

h 2 + 1 

h + 2 
- n 

dz 2 

1 3 ‘ d, 3 . 


z=0 


z=0 


h 3 + 
z=0 


(5) 


dE x (z) 


dz 


= jk cose 


Z=0 


d"E x (z) 


dz' 


- k 2 ( e - sin 2 © ) 
c 


z=0 


( 6 ) 


d J E x (z) 
dz 3 


z=0 


2 r . .2 de c 

k [ Dlc ( c c~ sin 8 ) cos8 + dz~ 


and c is the complex relative permittivity at the point z= 0 and 
dc c 

^ is the derivative of the permittivity with respect to z at the 
point z=0 + . Substituting equation (6) into equation (5) yields 


E^(h) = 1 - ( c + - sin 2 © ) - | k 2 h 3 e + ' 

2 2 3 3 

+ j kh cos© + k 2 —~ e + tan5 + - k 6 h ( e + - sin 2 ©) cos©j (7) 
where e is the real relative permittivity, e ' is the derivative 

T T 

of the real relative permittivity with respect to z, and tan5 + is 
the loss tangent. The subscript + denotes evaluation at z=0 + . 


The Taylor series expansion of the electric field about the 
point z=nh is written as 
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II II dE X ^ z ^ 

E x < z > = E x < nh > + dz^ 


( z - nh) + 
z=nh 


, 2^11 , , 

1 d E x < z > 

2 dz 2 


( z -nh) + • • • (8) 
z=nh 


Evaluating equation (8) at z=nh+h and z=nh-h and adding the 
result, thus giving 


II II 2 d E x ^ z ^ 

(nh+h) a 2E^ (nh) + h 


x 


dz' 


z=nh 


h 4 ^(z) 
12 dz 4 


- E^ I (nh-h) (9) 


z=nh 


where 


with 


d 2 E^(z) 

dz 2 


= k 2 ( c r (z) - sin 2 © ) E* J (z) 


z=nh 


z=nh 


d 4 E^(z) 


dz 


k 4 (e^z)-sin 2 © )E^ I (z) 


z=nh 


- 2k e r '(z) E* r (z) - kV'(z) E^ z (z) 


( 10 ) 


( 11 ) 


z = nh 


II' 

E* (nh) 


E x (nh) 


- E^(nh-h) 
h 


( 12 ) 


Employing the notation of Richmond, equation (9) , with equation 
(10) substituted, the electric field in the layer becomes 


Vl “ [ 2 - k2h2 < c n - sin2fl ) ] E „ - E n-1 + » 2 T§- (“) 

where E^ denotes equation (11) . The dielectric or plasma layer is 

divided into N segment so that Nh equals D, the thickness of the 
layer. The criteria for choosing the appropriate step size h is 
discussed in reference 7. The solution of equation (13) is 

determined in a step-by-step procedure as n varies from unity to 
N; i.e., for n=l the field at z=2h is computed from equation (13) 
since the fields at z=0 and z=h are given by equations (6) and (7) 
and the fourth derivative of the field at z=h is given by 

equations (10) -(12). By proceeding through the layer in this 
manner, the field at the region I and region II interface is 
determined. The field at this interface then enables one to 
determine the reflection coefficient. 
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Runcre-Kutta Method 


The equations needed for solving second-order differential 
equations of the form 


& = f(z,y) (14) 

dz 2 


by the Runge-Kutta method are given in reference 6. 
Equation (4) is written in this form as 


^ E (z) _ 2 II 

— 2 ~ k 2 ( e r (z) - sin 2 e )E xx (z) (15) 

dz 

Using the notation adapted earlier, the electric field at a 
general point z=nh where h is the spacing within the layer is 
represented as E n and at z=nh+h as E n+1 . From the fourth-order 

Runge-Kutta method, the solution takes the form 


E n+1 * E n + h 


[ ® ( 


k l + 2k 2 


)] 


E n+1 E n + 6 k l + 3 k 2 + 6 k 3 


(16) 

(17) 


where the primes denote the derivative with respect to z and 
k l - h £ ( nh ' E n) 

k 2 - h f( " h + i h, E n+ | E n ' + | kj 
k 3 = h f( nh + h, E n+ h E n ' + | k 2 ) 
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with 


f^nh, E n j = - k 2 ^ e(nh) - sin 2 © |E n 1 

f(nh+ | h, E n+ | e' + | kj =■ 

- k 2 [e£nh + §) - sin 2 e] [ V 5 E n + 5 k i] ' d' 

f(" h + h, E n+ h E n + | k 2 ) = 

- k 2 [e r (nh+ h) - sin 2 e] [ En + h E n '+ § kj 

J 

In this procedure, n varies from zero to N-l. The initial values 
/ 

E q and Eg are given by the first two equations in equation (6) . 
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DISCUSSION OF RESULTS 


Linear and exponential profiles were chosen to demonstrate 
the Richmond and Runge-Kutta methods for solving the differential 
equations. A comparison of these methods is made. These two 
profiles were selected since they also have exact solutions as 
shown in the Appendix. 

The particular linear profile chosen is shown in figure 2. 
The relative permittivity begins at a value of minus one (the 

exit point) and rises to a value of two at the front interface of 
the inhomogeneous layer. The real and imaginary parts of the 
electric fields at front interface for the Richmond (equation 
(13)) and the o Runge-Kutta (equation (16)) methods are shown in 
Table I for 5° increments in the angle of incidence. Also shown 
are the complex values of the electric fields for the exact 
solutions (Airy-equation (A-8)). The approximate computational 
time for each value in the Richmond method was 8 seconds compared 
to 12 seconds for the Runge-Kutta method. Both methods compare 
very favorably with the exact solutions as can be readily seen in 
Table I. 


TABLE I 


ELECTRIC FIELD DISTRIBUTION AT INTERFACE OF INHOMOGENEOUS 

LAYER-LINEAR PROFILE 


m= 3k and b=-l 



Richmond 

Runge- 

-Kutta 

Exact-Airy 

Angle 

Real 

Imag 

Real 

Imag 

Real 

Imag 

0 

0.9918 

0.9178 

0.9919 

0.9179 

0.9919 

0.9179 

5 

0.9954 

0.9155 

0.9955 

0.9156 

0.9955 

0.9156 

10 

1.0063 

0.9086 

1.0065 

0.9086 

1.0065 

0.9086 

15 

1.0243 

0.8969 

1.0244 

0.8969 

1.0244 

0.8969 

20 

1.0488 

0.8800 

1.0489 

0.8800 

1.0489 

0.8800 

25 

1.0794 

0.8577 

1.0794 

0.8577 

1.0794 

0.8577 

30 

1.1150 

0.8296 

1.1152 

0.8296 

1.1152 

0.8296 

35 

1.1550 

0.7953 

1.1552 

0.7953 

1.1552 

0.7953 

40 

1.1983 

0.7543 

1.1985 

0.7543 

1.1985 

0.7543 

45 

1.2435 

0.7065 

1.2437 

0.7065 

1.2437 

0.7065 

50 

1.2894 

0.6515 

1.2895 

0.6516 

1.2895 

0.6516 

55 

1.3345 

0.5895 

1.3336 

0.5895 

1.3346 

0.5895 

60 

1.3775 

0.5206 

1.3775 

0.5206 

1.3775 

0.5206 

65 

1.4166 

0.4452 

1.4167 

0.4452 

1.4167 

0.4452 

70 

1.4508 

0.3640 

1.4509 

0.3640 

1.4509 

0.3640 

75 

1.4788 

0.2777 

1.4789 

0.2777 

1.4789 

0.2777 

80 

1.4995 

0.1874 

1.4997 

0.1874 

1.4997 

0.1874 

85 

1.5124 

0.0944 

1.5125 

0.0944 

1.5125 

0.0944 

90 

1.5167 

0.0000 

1.5168 

0.0000 

1.5168 

0.0000 
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The exponential profile selected for comparison is shown in 
figure 3. The relative permittivity e r begins at a value of one 

and rises to a value of e at the front interface of the 
inhomogeneous layer. The values of the real and imaginary parts 
of the electric fields at the front interface for the Richmond and 
Runge-Kutta methods are shown in Table II for 5° increments in the 
angle of incidence. Also shown are the complex values of the 
electric fields for two angles of incidence ( 0°and 30°) for the 
exact solution (Bessel-equation (A-17 )); these two angles were 
selected because, for the parameters used in these computations, 
they yield integer order Bessel functions which are easily 
calculable. The approximate computational time for each value in 
the Richmond method was 9 seconds compared to 12 seconds for the 
Runge-Kutta method. Both methods compare very favorably with each 
other as well as to the two exact solutions as can be readily seen 
in Table II. 


TABLE II 


ELECTRIC FIELD DISTRIBUTION AT INTERFACE OF INHOMOGENEOUS 

LAYER-EXPONENTIAL PROFILE 


§=k 



Richmond 

Runge- 

-Kutta 

Exact- 

-Bessel 

Angle 

Real 

Imag 

Real 

Imag 

Real 

Imag 

0 

0.3736 

0.7411 

0.3738 

0.7411 

0.3738 

0.7411 

5 

0.3766 

0.7393 

0.3767 

0.7393 



10 

0.3851 

0.7340 

0.3853 

0.7340 



15 

0.3993 

0.7249 

0.3994 

0.7249 



20 

0.4187 

0.7119 

0.4188 

0.7119 



25 

0.4427 

0.6945 

0.4430 

0.6945 



30 

0.4711 

0.6725 

0.4713 

0.6725 

0.4722 

0.6722 

35 

0.5028 

0.6455 

0.5030 

0.6455 



40 

0.5372 

0.6130 

0.5373 

0.6130 



45 

0.5731 

0.5748 

0.5733 

0.5749 



50 

0.6096 

0.5308 

0.6098 

0.5308 



55 

0.6456 

0.4809 

0.6457 

0.4809 



60 

0.6798 

0.4251 

0.6800 

0.4252 



65 

0.7112 

0.3639 

0.7113 

0.3640 



70 

0.7385 

0.2976 

0.7387 

0.2978 



75 

0.7610 

0.2273 

0.7611 

0.2273 



80 

0.7776 

0.1535 

0.7777 

0.1535 



85 

0.7878 

0.0774 

0.7880 

0.0774 



90 

0.7912 

0.0000 

0.7914 

0.0000 




Reflection coefficients at the front interface of the 
inhomogeneous layers are also of interest. The equations required 
for computing the reflection coefficients are given in reference 
7. Magnitude and phase of the reflection coefficients for the 
linear and exponential profiles discussed in this paper are shown 
in Tables III and IV. 
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TABLE III 


REFLECTION COEFFICIENTS AT INTERFACE OF INHOMOGENEOUS 
LAYER-LINEAR PROFILE 


m= 3k and b=-l 



Richmond 

Runge-Kutta 

Angle 

Magnitude 

Phase 

in radians 

Magnitude 

Phase 

in radians 

0 

0.2925 

-.0548 

0.2929 

-.0825 

5 

0.2932 

-.0480 

0.2936 

-.0757 

10 

0.2954 

-.0277 

0.2957 

-.0553 

15 

0.2991 

0.0030 

0.2992 

-.0211 

20 

0.3044 

0.0542 

0.3043 

0.0268 

25 

0.3116 

0.1159 

0.3112 

0.0889 

30 

0.3208 

0.1919 

0.3201 

0.1655 

35 

0.3324 

0.2828 

0.3314 

0.2571 

40 

0.3473 

0.3897 

0.3457 

0.3644 

45 

0.3659 


0.3638 

0.4885 

50 

0.3896 


0.3870 

0.6308 

55 

0.4202 


0.4169 

0.7934 

60 

0.4604 


0.4564 

0.9797 

65 

0.5137 


0.5092 

1.1946 

70 

0.5862 


0.5811 

1.4463 

75 

0.6840 


0.6788 

1.7485 

80 

0.8090 


0.8048 

2.1221 

85 

0.9384 


0.9367 

2.5901 

90 

1.0000 


1.0000 

3.1416 


10 












TABLE IV 


REFLECTION COEFFICIENTS AT INTERFACE OF INHOMOGENEOUS 
LAYER-EXPONENTIAL PROFILE 


j3=k 



Richmond 

Runge-Kutta 

Angle 

Magnitude 

Phase 

in radians 

Magnitude 

Phase 

in radians 

0 

0.2788 

-2.4876 

0.2825 

-2.4801 

5 

0.2800 

-2.4863 

0.2837 

-2.4791 

10 

0.2836 

-2.4833 

0.2874 

-2.4762 

15 

0.2898 

-2.4780 

0.2935 

-2.4714 

20 

0.2986 

-2.4710 

0.3024 

-2.4651 

25 

0.3105 

-2.4631 

0.3144 

-2.4575 

30 

0.3258 

-2.4540 

0.3297 

-2.4494 

35 

0.3449 

-2.4450 

0.3489 

-2.4414 

40 

0.3686 

-2.4372 

0.3727 

-2.4345 

45 

0.3979 

-2.4320 

0.4021 

-2.4302 

50 

0.4340 

-2.4311 

0.4382 

-2.4303 

55 

0.4784 

-2.4373 

0.4828 

-2.4374 

60 

0.5334 

-2.4542 

0.5376 

-2.4553 

65 

0.6009 

-2.4868 

0.6051 

-2.4889 

70 

0.6831 

-2.5427 

0.6969 

- 2 . 5455 

75 

0.7791 

-2.6306 

0.7823 

-2.6338 

80 

0.8808 

-2.7601 

0.8828 

-2.7631 

85 

0.9656 

-2.9348 

0.9663 

-2.9367 

90 

1.0000 

-3.1416 

1.0000 

-3 . 1416 
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CONCLUDING REMARKS 


The numerical techniques of Richmond and Runge-Kutta were 
shown to be quite efficient for solving the second-order 
differential equation which describes the field within a plane 
inhomogeneous layer of permittivity, upon which a perpendicularly- 
polarized plane wave has been assumed incident. Richmond's method 
seems to be little more efficient than the Runge-Kutta method, at 
least for the two relative permittivity profiles considered. 

Although the techniques of Richmond and Runge-Kutta were 
applied to analytical profiles, they should be equally valid for 
arbitrary relative permittivity profiles or even profiles which 
are known discretely. This is important because realistic 
relative permittivity profiles may take on almost any shape. 
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APPENDIX 


Exact solutions of equation (4) are available for only a few 
relative permittivity profiles [3,4]. A brief derivation of the 
solutions of equation (4) for linear and exponential profiles is 
given here. 

For the linear profile, the relative permittivity is assumed 
to vary as 


e r (z) = mz + b 


(A-l) 


where m is the slope of the profile and b is the beginning value 
of the profile ( at z=0 ) . Substituting equation (A-l) into 

equation (4) yields 

d2 V z > 2 2 

= + k ( mz + b - sin e ) E v (z) = 0 (A-2) 

dz^ x 


By making the substitution 

z = Mv + B (A- 3 ) 

equation (A-2) becomes 


d"E x (v) 

dv 2 


| M 2 k 2 £ mMv + mB + b- sin 2 © Jj E x (v) = 


By equating the term in braces to -v and solving for 
is found that 


0 (A-4 ) 

M and B, it 


d 2 E x (v) 

2 V M V ) = 0 (A-5) 

dv x 


where 


1 



. 2 

_ sin 9 - b 


From equations (A-3) and (A-6), it follows that 


(A-6 ) 
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v = - ( mz + b- sin 2 © ) | jjj j 


(A-7) 


The solutions of equation (A-5) are the two independent Airy 
functions, Ai(y) and Bi(v) [4]. The general solution of equation 
(A-5) , therefore is written as 


E x (v) = C^Ai (v) + C 2 Bi(v) 


(A-8) 


where C. and C_ are constants. These constants are determined by 

applying the appropriate boundary conditions at z=0. The boundary 
conditions are given by the first two equations in equation (6) 
but which must be expressed in terms of the new independent v. 
The evaluation of the constants is given in matrix form as 


f C, 


1 

0 


Bi(vo) - Bi(vo) ) 
/ 

- Ai(vo) Ai(vo) 


-j | ^ ) 3 cos© 
1 m 


( A-9 ) 


where 


0= Ai(vo)Bi(Vo) - Ai(vo)Bi(Vo) 


Vo = - 




k 

m 


. 2 

b - sin 9 


(A-10 ) 


and the primes denote the derivative with respect to v. For the 
solutions in this paper m=3k and b=-l. Once the C's are determined, 
the electric field at any point in the layer is computed by using 
equation (8) . 

For the exponential profile, the relative permittivity is 
assumed to vary as 


e r (z) = e^ z 


(A-ll) 


Substituting equation (A-ll) into equation (4) yields 


d E *( z > 2 
— ^ + k 2 

dz 


e* 32 - sin 2 © 1 E (z) = 0 
; X 


(A-12 ) 


Following the substitution given in reference 4, 
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0Z 

V = 2k e 2 

0 ® 

the differential equation is now written as 

, 2 , 


2 d E x (v) dE x (v) 


+ v 


dv‘ 


dv 


+ ( V 2 - p 2 ) E x (v) = 0 


(A-13) 


(A-14) 


where 


p 2 = ikl s . n 2 e 

0 


(A-15 ) 


If p is zero or a positive integer, the general solution of 
equation (A-14) becomes 


V V > “ B l J p< v > + B 2 Y p< v > 


(A— 16) 


where Jp(v) is the Bessel function of the first kind, Y^fv) is the 

Bessel function of the second kind, and the B's are constants. 
Upon applying the appropriate boundary conditions at z=0, the 
matrix for the B's becomes 


(A— 17 ) 



1 

" Y P < V “) 



r 1 

■ B 2 - 

D 

- J p (V.) 

J p (Vo) 


j cose 

^ j 


where 


D= J (Vo)YJVo) - J (Vo) Y (Vo) 


Vo = 


P 

2k 

0 


P 


(A-18 ) 


and the primes denote the derivative with respect to v. For the 
solutions in this paper, /3=k with kD=l. Once the B's are 
determined, the electric field at any point in the layer is 
computed by using equation (16) . If p is not zero or a positive 
integer, the general solution of equation (A-14) becomes 


V V > - D l J p< V > + D 2 J -p< V > 


(A— 19 ) 


where the constants D 1 and D 2 are determined in the manner just 
described. 
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Figure 1 .-Perpendicularly-polarized plane wave incident 
on an inhomogeneous plasma layer. 
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Figure 3. -Exponential profile. 
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